function [stayf0,staym0,stayf1,staym1] = stayM(H, wf, wm, t, det_earn_f, det_earn_m, home_prod, psi, th1, th_tau, lambda, dt)

W0=w(det_earn_m,t,H-1,wm,dt)+w(det_earn_f,t-(H-1),0,wf,dt); 
W1=w(det_earn_m,t,H-1,wm,dt)+psi*w(det_earn_f,0,0,0,dt); 

stayf0=log(lambda*(W0/2)^2) + th1 + th_tau; 
staym0=log((1-lambda)*(W0/2)^2) + th1 + th_tau;
stayf1=log(lambda*(W1/2)^2) + th1 + th_tau; 
staym1=log((1-lambda)*(W1/2)^2) + th1 + th_tau;

end
